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ABSTRACT 

The excursion set approach allows one to estimate the abundance and spatial dis- 
tribution of virialized dark matter haloes efficiently and accurately. The predictions 
of this approach depend on how the nonlinear processes of collapse and virialization 
are modelled. We present simple analytic approximations which allow us to compare 
the excursion set predictions associated with spherical and ellipsoidal collapse. In par- 
ticular, we present formulae for the universal unconditional mass function of bound 
objects and the conditional mass function which describes the mass function of the 
progenitors of haloes in a given mass range today. We show that the ellipsoidal col- 
lapse based moving barrier model provides a better description of what we measure 
in the numerical simulations than the spherical collapse based constant barrier model, 
although the agreement between model and simulations is better at large lookback 
times. Our results for the conditional mass function can be used to compute accurate 
approximations to the local-density mass function which quantifies the tendency for 
massive haloes to populate denser regions than less massive haloes. This happens be- 
cause low-density regions can be thought of as being collapsed haloes viewed at large 
lookback times, whereas high-density regions are collapsed haloes viewed at small 
lookback times. Although we have only applied our analytic formulae to two simple 
barrier shapes, we show that they are, in fact, accurate for a wide variety of moving 
barriers. We suggest how they can be used to study the case in which the initial dark 
matter distribution is not completely cold. 
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1 INTRODUCTION 

Bond et al. (1991) described an approach which allowed 
them to use the statistical properties of the initial density 
fluctuation field to derive an estimate of the number density 
of collapsed dark matter haloes at later times: the so-called 
universal, unconditional mass function (Press & Schechter 
1974). Lacey & Cole (1993) showed how this model could 
be extended to estimate the rate at which small objects 
merge with each other to produce larger objects. This al- 
lowed them to estimate the conditional mass function of 
subhaloes within parent haloes. Lacey & Cole also provided 
formulae for the distribution of halo formation times; Nusser 
& Sheth (1999) provided formulae for the distribution of the 
halo mass at formation. Sheth (1996) and Sheth & Pitman 
(1997) showed how various higher order statistical proper- 
ties of the forest of merger history trees associated with the 
formation of these objects could also be estimated within 
this approach. Mo & White (1996) and Sheth & Lemson 



(1999a) showed how information about the forest of merger 
history trees could be used to quantify the extent to which 
dark haloes are biased tracers of the matter distribution. Mo, 
Jing & White (1997) provided predictions for the higher or- 
der moments of the spatial distribution of the haloes, and 
Sheth (1998) showed how to use the approach to estimate 
the probability that a randomly placed cell contains a cer- 
tain amount of mass. Clearly, the approach has been very 
useful. 

The approach combines the simple physics of the spher- 
ical collapse model with the assumption that the initial 
fluctuations were Gaussian and small. The problem of es- 
timating any one of the quantities listed above is reduced 
to solving a problem associated with the crossing of an ap- 
propriately chosen barrier by particles undergoing Brownian 
motion; the Brownian nature of the motion comes from the 
assumption of Gaussian initial conditions, and the barrier 
shape is specified by the spherical collapse model (e.g. Sheth 
1998). Hence, this is often called the excursion set approach. 
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Tormen (1998) reported that the spherical collapse 
based excursion set predictions did not describe the con- 
ditional mass function of subclumps in his simulations well. 
Motivated by this, Sheth, Mo & Tormen (2001) discussed 
a simple way of modifying the excursion set approach to 
incorporate the effects of ellipsoidal, rather than spherical 
collapse. On average, initially denser regions collapse before 
less dense ones. This means that, at any given epoch, there 
is a critical density which must be exceeded if collapse is to 
occur. In the spherical collapse model, this critical density 
does not depend on the mass of the collapsed object. How- 
ever, in their parametrization of ellipsoidal collapse, Sheth, 
Mo & Tormen showed that, of the set of objects which col- 
lapse at the same time, the less massive ones must initially 
have been denser than the more massive ones, since the less 
massive ones would have had to hold themselves together 
against stronger tidal forces. They argued that this could 
be incorporated into the excursion set approach, simply and 
effectively, if not rigorously, by changing the barrier shape. 
In essence, whereas the barrier associated with spherical col- 
lapse is one whose height does not depend on distance from 
the origin of the walk, the one associated with ellipsoidal dy- 
namics increases with distance. They showed that the excur- 
sion set approach with a moving barrier was able to provide 
a good fit to the universal halo mass function. 

This paper is devoted to a more detailed discussion of 
moving barrier models. In general, moving barrier models 
have a richer structure than the constant barrier model. 
For example, the approach with spherical dynamics predicts 
that, at any given time, all the mass in the universe is bound 
up in collapsed objects, whereas a small fraction of the mass 
remains unbound in the case of ellipsoidal dynamics. In ad- 
dition, whereas clustering is strictly hierarchical in the case 
of spherical dynamics, incorporating ellipsoidal collapse into 
the excursion set approach results in a model in which frag- 
mentation as well as mergers may occur — the approach pre- 
dicts that some small haloes fragment before they are sub- 
sumed into larger ones. Appendix A describes some of these 
features, which may (or may not!) provide better approxi- 
mations to the physics of graviational instability than does 
the constant barrier model, in more detail. The Appendix 
also describes how the results of this paper can be used to 
model the halo mass function in warm dark matter scenarios 
such as that revisited by Bode, Ostriker & Turok (2001). 

Section ^ shows how moving barrier models can be used 
to make simple analytic estimates of a number of statistical 
quantities which are routinely measured in numerical cos- 
mological simulations. The primary results of Section ^ are 
equations (^) and (^), which are accurate for a large class 
of moving barrier shapes. To illustrate how these formulae 
work, we show the result of inserting the ellipsoidal collapse 
moving barrier of Sheth, Mo & Tormen (2001) into these 
formulae. Section 



2.1 



presents the number density of bound 
objects as a functi on o f mass (the unconditional mass func- 
tion), and Section [2!^ describes a simple efficient algorithm 
for generating it. Section |2.3| presents the average number 
of progenitor subhaloes as a function of subhalo mass, for a 
wide range of specified parent halo masses (the conditional 
mass function), and Section 2A shows how the halo mass 



ing barrier shape, is quite accurate. Section |2.5| shows that, 
at small lookback times, neither the constant nor the moving 
barrier predictions describe the conditional mass functions 
in the simulations particularly well, although the agreement 
at large lookback times is quite good if the spherical collapse 
constant barrier is used, and even better if the ellipsoidal 
collapse moving barrier is used. 

Section ^ shows the result of considering more compli- 
cated moving barrier excursion set models. In particular, it 
shows the result of considering the full six-dimensional ran- 
dom walk associated with Gaussian random fields, rather 
than the one-dimensional simplification proposed by Sheth, 
Mo & Tormen (2001). It shows that their simplification is 
actually quite accurate. Details associated with the calcula- 
tions in this section are presented in Appendix ^. 

Section 1^ discusses some simple implications of our find- 
ings. Although, in this paper, we concentrate on the moving 
barrier derived by Sheth, Mo & Tormen (2001), we think 
it worth stressing that our analytic formulae are more gen- 
eral: they are accurate for a wide variety of moving barrier 
shapes. 



2 THE MOVING BARRIER MODEL 

As discussed in the introduction, we will mainly be inter- 
ested in the first crossing distributions of uncorrelated Brow- 
nian motion random walks. Following Bond et al. (1991) 
and Lacey & Cole (1993), these first crossing distributions 
can be used to provide useful approximations to what have 
come to be called the conditional and unconditional mass 
functions of the dark halo distribution. The results of this 
section should be thought of as generalizations of the re- 
sults in Lacey & Cole (1993). Whereas they restricted their 
attention to a barrier of fixed height, this section presents 
analytic formulae which approximate the barrier crossing 
distribution for a wide class of moving barriers. 

Before we begin, we think a word on notation is helpful. 
We have chosen to present our formulae for the first crossing 
distributions using the same notation as Lacey and Cole 
(1993). This means that we use the symbol S to represent 
the variance in the density fluctuation field when smoothed 
on some scale, which is usually denoted . However, some 
of our formulae can be written in terms of the scaled variable 
(cr/(7,)^ = S/St, for some suitably defined St. When this is 
possible, we also write our formulae in terms of 1/ = S, / S . 



2.1 The unconditional mass function 

To illustrate how our formulae work, rather than use the 
spherical collapse barrier of constant height, we will use 
the moving barrier shape derived by Sheth, Mo & Tormen 
(2001) — the one associated with ellipsoidal collapse: 



B(<7^ z) = ^^(5,c(^) [l + l3 (ai^)-" 



(1) 



function depends on the surrounding density field (the local- 
density mass function). Comparison with simulations shows 
that the approach, with the ellipsoidal collapse based mov- 



where ly = [5Bc{z)/o-{m)]^ . In the expression above, Sbc{z) is 
the critical overdensity required for spherical collapse at z, 
extrapolated using linear theory to the present time (it is 
approximately 1.686(1 + z) if fl = 1), and cr(m) is the rms 
of the initial density fluctuation field when it is smoothed 
on a scale which contains mass m, extrapolated using linear 
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theory to the present time. The parameters /3 ~ 0.485 and 
a « 0.615 come from ellipsoidal dynamics (the spherical 
collapse model has a = and /3 = 0), and the value a ~ 0.7 
comes from normalizing the model to simulations (as we 
discuss later, it may be more accurate to set a ~ 0.75 or so). 
The spherical collapse model sets B{a^, z) = 5sc{z)\ because 
B is then independent of a, spherical collapse is said to be 
associated with a barrier of constant height. 

As Sheth, Mo & Tormen noted, B{a^,z) scales simi- 
larly to the natural scaling associated with random walks: 
multiplying the barrier height by c and rescaling a by the 
same factor results in a barrier of the same shape. They 
exploited this property as follows: first, they studied what 
happens when z — 0. Whereas an analytic expression for 
the first crossing distribution of a barrier of constant height 
has been known for some one hundred years or so, there 
is as yet no analogous solution for the ellipsoidal collapse 
moving barrier. (One of the main results of this section is to 
provide a good analytic approximation to this solution.) So, 
Sheth, Mo & Tormen simulated a large number of random 
walks and computed the distribution of first crossings, f{S), 
of the barrier B{a^ — S, z = 0) numerically. The scaling of 
the barrier shape meant that their numerical solution could 
be scaled to represent the first crossing distribution at any 
other time also. Therefore, when providing an analytic fit to 
their simulated first crossing distribution, they expressed it 
in scaled variables: 



ufiu)=A 



1 + {aiyy 



av\ ^1"^ e 
~2 



-au/2 



(2) 



where / is the distribution of first crossings, u and a axe 
the variables defined above, p — 0.3 and A is determined by 
requiring that the integral of /{v) over all u give unity. The 
distribution associated with the spherical collapse constant 
barrier is got by setting a = 1, p = and A — 1/2. 

In the excursion set approach, the average comoving 
number density of halos of mass m, often called the universal 
or unconditional halo mass function n{m,z), is related to 
the first crossing distribution (this is why the shape of the 
barrier influences the shape of the mass function) by 



2 n(m, z) d In m 
p d In 1/ ' 



(3) 



where p denotes the average comoving density of the back- 
ground. Because the first crossing distribution evolves self- 
similarly, so does the halo mass function. Fig. 2 in Sheth & 
Tormen (2001) shows that, in the GIF simulations of clus- 
tering in SCDM, OCDM and ACDM cosmologies, n(m, z) is 
well fit by this expression (at least over the range z = to 
2 = 4). That is to say, the excursion set prediction that the 
universal unconditional mass function evolves self-similarly, 
is actually a very good approximation. 

Although this fitting function (equation ^) for the first 
crossing distribution is extremely useful, when we consider 
conditional mass functions in the next section, it will turn 
out that we need to compute a new fitting function for each 
parent mass range of interest. This is, in principle, a compu- 
tational bottleneck. For this reason, we think it more useful 
to provide a formula for the first crossing distribution which 
is more amenable for use in what is to follow. 

We have analytic formulae for the first crossing distri- 
bution only in the case of constant (e.g. Bond et al. 1991) 
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Figure 1. First crossing distributions and the universal uncondi- 
tional halo mass function. Histogram shows the distribution ob- 
tained by simulating random walks which are absorbed on the 
ellipsoidal collapse moving barrier; solid curve shows our analytic 
approximation to this distribution (equation Dashed curve 
shows the distribution which fits the halo mass function in cos- 
mological simulations (equation and dotted curve shows the 
distribution associated with spherical collapse. 



and linear (Sheth 1998 and Appendix A of this paper) bar- 
riers. These solutions suggest the following approximation 
which we have found to work rather well. For a wide range 
of moving barrier shapes B{S), the first crossing distribution 
is well approximated by 

f{S) dS = \TiS)\ exp (-^) (4) 



where T{S) denotes the sum of the first few terms in the 
Taylor series expansion of B(S): 

n = 

Notice that this expression gives the correct answer in the 
case of constant and linear barriers (in which only the first, 
or the first two terms of the series are non-zero). For the 
ellipsoidal barrier shape, we find we get reasonable accu- 
racy to the numerical result if we truncate this expansion 
at n = 5. The accuracy of this formula increases as the dis- 
tance between the start of the walk and the barrier height 
at that initial position, increases (this distance is ^/a for the 
ellipsoidal collapse barrier of equation |^) . 

To illustrate. Fig. ^ shows the result of setting the bar- 
rier B{S) to equal that given by equation (|l|) at z = 0, 
then generating a large ensemble of random walks, and so 
constructing the first crossing distribution. The histogram 
shows this Monte-Carlo distribution, and the solid curve 
shows the approximation presented in equation (W): 



"^'^g-ai/[l + /3(ai^) 

27r" 



1 + 
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a(a — 1) 

- + ■ 



/2 



(5) 



These should be compared with the distribution of first 
crossings of a barrier of constant height set equal to Ssc {z = 
0) (dotted curve), and equation (^) which fits the GIF sim- 
ulations (dashed curve). With the exception of the dotted 
curve (the one associated with the constant barrier spheri- 
cal collapse model) the other three curves are in reasonably 
good agreement. We will exploit this fact in the next subsec- 
tion, when we derive a simple expression for the conditional 
mass function associated with ellipsoidal collapse. 

Before we do so, however, we think it useful to point 
out that there are some generic properties of moving barrier 
models which arise from the scaling of the barrier shape. The 
fact that the barrier shape scales in the way it does means 
that moving barrier models in which the barrier height in- 
creases with S are somewhat more complicated than the 
constant barrier model, so they may be used to model a 
wider variety of physical processes. For example, our final 
expression for the moving barrier crossing distribution is not 
normalized to unity (notice that the solid curve is always be- 
low the dashed one). This is because, if the barrier height 
increases more steeply than linearly with a then not all walks 
cross the barrier. Appendix A uses a simple analytic moving 
barrier model to illustrate some of these features. 



2.2 Generating the unconditional mass function 

Before we move on to study the conditional mass function, 
we thought it useful to describe a fast and simple algorithm 
for generating random numbers which can be used to pro- 
vide the correct distribution of halo masses. To generate 
numbers drawn from the spherical collapse mass function 
(equation ^ with A = 1/2 and p = 0) is easy because, in this 
case, f{u) can be got by generating a Chi-squared random 
variate. This means that, in the case of spherical collapse, 
the unconditional mass function can be generated quickly 
and easily because generating Gaussian random variates is 
easy. In particular, one simply generates a Gaussian random 
variate x, and then sets au — . 

The ellipsoidal collapse mass function (equation ^ with 
A = 0.32218 and p = 0.3) cannot be transformed to a Gaus- 
sian, so constructing an algorithm for generating it is not 
so straightforward. We have found that first generating a 
Gaussian variate x, and then setting av = |a;|"^ ''/(l -I- la::!"'^ '') 
is accurate to within a percent or so over the range 0.01 < 
u < 100. The speed of this algorithm compensates for the 
fact that it does not exactly produce variates drawn from 
the ellipsoidal collapse mass function. 



2.3 The conditional mass function 

As stated above, we do not know of an analytic expres- 
sion for the first crossing distribution associated with bar- 
riers which have the form given in equation (|l]). How- 
ever, we do have two reasonably accurate fitting formulae — 
equations (^) and (^ — to this distribution. One might have 
thought that we could use them to make an estimate of the 
conditional mass function as follows. 
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Figure 2. Examples of random walks used to construct the con- 
ditional mass functions associated with the ellipsoidal collapse 
moving barriers at z = (dotted curve) and 2 = 2 (solid curve). 



Bond et al. (1991) and Lacey & Cole (1993) argued that 
conditional mass functions could be estimated by consider- 
ing the successive crossings of boundaries associated with 
different redshifts. The first crossing of two constant barri- 
ers of different heights has an analytic solution, so they were 
able to provide analytic estimates for the conditional mass 
function associated with the spherical collapse model. Such 
a formula is very useful, because, once the conditional mass 
function is known, the forest of merger history trees can be 
constructed using the algorithm described by Sheth & Lem- 
son (1999b), from which the nonlinear stochastic biasing as- 
sociated with this mass function can be derived using the 
logic of Mo & White (1996) and Sheth & Lemson (1999a). 

In the constant barrier case, the conditional mass func- 
tion is computed by considering walks which start from 
[cr^ — S,5bc{zo)] rather than from the origin, and then in- 
tersect the constant barrier Ssc{zi) at, say, s. This is eas- 
ily computed, because, despite the shift in the origin, the 
second barrier is still one of constant height. Since the 
first crossing distribution of a barrier of constant height is 
known, (recall it was just such a distribution which was as- 
sociated with the universal, unconditional mass function), 
the conditional mass function can also be written analyti- 
cally. Essentially, f{s\S) has the same form as the uncondi- 
tional mass function f{S), but with the change of variables: 
5sc{z) Ssc{zi) — Sac{zo) and S ^ s — S. 

One might have wondered if the same change of vari- 
ables in equation (^ provides a good description of the con- 
ditional mass function associated with the ellipsoidal col- 
lapse moving barrier. Unfortunately, because the barrier 
shape is not linear in S, changing the origin of the coor- 
dinate system does not yield a barrier of exactly the same 
functional form. Specifically, the shape of 

B{s,zi)- B{S,zo) = VE5i[l + Ps°'/{aSl) 



where 5i = 5{zi) etc., can be written as a constant plus 
a term which scales as (s — S)" only when a equals zero 
or one. This means that, formally, the solution to the two 
barrier problem associated with ellipsoidal dynamics is not 
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given by a simple rescaling of the unconditional ellipsoidal 
collapse mass function. Therefore, we cannot simply rescale 
the fitting function of equation (|^) to get a reliable estimate 
of the conditional mass function: the two barrier problem 
associated with moving barriers must, in general, be solved 
numerically. 

This is discouraging because it means that, in principle, 
we must find a different fitting function for each choice of 
condition, because each condition corresponds to a differ- 
ent origin, say, {Bo, So), and so to a slightly different bar- 
rier shape. Of course, the result can be generated relatively 
quickly in at least two ways. The first is to solve the integral 
equation associated with this barrier numerically (Monaco 
1997b; Sheth 1998). The second is to simply simulate the 
random walk trajectories and so construct the first crossing 
distribution directly. 

Fig. ^ shows an example of what is involved in solv- 
ing the two barrier problem numerically in this way. The 
smooth dotted and solid curves show the ellipsoidal col- 
lapse barrier B{S,z) of equation (Q) scaled to 2; = and 
z = 2, respectively. Jagged curves show a few representative 
random walk trajectories: they start at the barrier position 
B{S, 2 = 0), where S{M) is given by the GIF SCDM power 
spectrum, and M/M* = 2. These random walks are followed 
until they first cross the barrier B{s,z = 2). The value of 
5* at which this happens is stored and used to make plots 
like those shown below. The random walks were generated 
by the same Monte-Carlo code that was used to generate 
Fig. except that there the barrier shape was given by 
equation ( |Al[ ). Fig. Al shows that this Monte-Carlo code 
works correctly. One possible approach to the excursion set 
conditional mass function, then, is to simply generate it nu- 
merically, as the need arises. 

Before providing a detailed comparison of the condi- 
tional mass functions generated using this Monte-Carlo 
model and those in numerical simulations, it is useful to 
study a simple limiting case. Fig. |2| shows that for S/St < 
0.5, the height of the barrier is approximately constant. At 
small S, the only difference between the barrier at two red- 
shifts, and the spherical collapse constant barriers, arises 
from the factor of a = 0.707. This has the following conse- 
quence. At small lookback times (small redshift differences), 
most random walk trajectories will intersect the barrier be- 
fore they have travelled very far along the 5* axis. For these 
trajectories, the barrier is, effectively one of constant height. 
This means that the conditional mass function for massive 
haloes at small lookback times will have the same shape as 
that predicted by the constant barrier, with one small dif- 
ference. The factor of a = 0.707 has the effect of slightly 
reducing (by a factor of ^/a) the effective redshift difference 
relative to the original constant barrier model. As a result, 
the GIF barrier suggests that massive haloes at small look- 
back times will be slightly more massive than the original 
constant barrier model predicts. Since \/0.707 is close to 
unity, this effect will be small. In practice, we only expect 
the barrier predictions to differ significantly from those of 
the constant barrier for small haloes, or at large lookback 
times (high redshift). This is encouraging, because these are 
precisely the regimes in which simulations suggest that the 
constant barrier model is inaccurate (Tormen 1998). 

In addition to Monte-Carloing the conditional mass 
functions, we can use the results of the previous subsec- 



tion to derive a simple analytic expression for their shape. 
We can do this because our barrier crossing formula is rea- 
sonably accurate for a rather wide range of barrier shapes. 
A glance at Fig. |^ shows that the barrier shapes associated 
with the conditional mass functions are not likely to be too 
different from those associated with the unconditional func- 
tion, so we should be able to use equation (^ to approximate 
most of the conditional mass functions we will be interested 
in. In practice, this can be done by simply making the ap- 
propriate replacements B — » B{s) — B{S) and S ^ s — S 
in equation (^. At the risk of being repetitive, our approx- 
imation for the conditional mass function is: 



J^(m,5i\M,5o) dm = — f{m, Si\M, So) dm 



(6) 



where f{m\M) dm = f{s\S) ds with 



fis\S) ds 



and 



\T{s\S)\ 



exp 



[B{s)-BiS)] 
2(5-5) 



ds 



(7) 



T{s\S)=Yl 

n=0 



(s- 



d"[B{s)-B{S)] 



Fig. ^ compares this approximation with the actual nu- 
merical Monte-Carlo distribution, and compares both with 
the actual distribution measured in the cosmological sim- 
ulations. (We will show a comparison with what one gets 
by rescaling equation ^ shortly.) The histograms (without 
error bars) show the conditional mass functions generated 
using our Monte-Carlo code, for parent haloes in the mass 
range 1 < M/M* < 2 (upper curves) and 16 < M/M, < 32 
(lower curves); the upper curves have been shifted upwards 
by a factor of ten for clarity. The smooth curves show the 
associated analytic approximation; they describe the results 
of our numerical walks reasonably well. Therefore, in the 
comparisons to follow, we will sometimes only show the an- 
alytic curves. Recall that the analytic formula should be 
most accurate for the high-redshift progenitors of massive 
parents, and least accurate when the parents are not very 
massive. The figure shows some evidence that this is true. 
Both these curves should be compared with the symbols 
which have error bars — they show the conditional mass func- 
tions measured in the ACDM simulations. The figure is quite 
encouraging: our excursion set predictions are in reasonably 
good agreement with the cosmological simulations. In addi- 
tion, the figure shows that the analytic approximation (equa- 
tion 1^ is reasonably accurate. 

Fig. 1^ shows a similar comparison in the case of SCDM. 
In this case, we have chosen to not show the random walk 
histogram, since it is quite well fit by our formula. Instead, 
the figure shows the conditional mass functions measured in 
the SCDM simulations (symbols with error bars) and the 
various analytic approximations discussed earlier (smooth 
curves). In order of accuracy, these are equation (^) (solid), 
the distribution associated with rescaling the unconditional 
mass function of equation (^) (dashed) , and the conditional 
distribution associated with the constant barrier, spherical 
collapse model of Lacey & Cole (1993) (dotted). (The finite 
mass resolution of the simulations means that the distribu- 
tions are artificially truncated at low masses.) 

Notice that the spherical collapse based dotted curves 
are often quite different from the N-body simulation sym- 
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Figure 3. Excursion set conditional mass functions at z, for parent haloes identified at the present time. The two sets of curves in 
each panel are for parent haloes in the mass range 1 < M/M, < 2 (upper curves) and 16 < M/M* < 32 (lower curves); the upper 
curves have been shifted upwards by a factor of ten for clarity. Symbols with error bars show the distributions measured in the ACDM 
simulations, histograms show the result of generating the first crossing distribution by simulating an ensemble of 10* random walks, 
smooth solid curves show the analytic approximation discussed in the text, and dotted curves show the distribution associated with 
barriers of constant height. 



bols. This discrepancy is similar to that first noticed by Tor- 
men (1998); haloes in the simulations seem to hold them- 
selves together at higher redshift than the spherical collapse 
model predicts. Notice also that, whereas the dashed curves 
we get by rescaling the unconditional halo mass function 
are certainly better fits to the cosmological simulations, the 
solid curves, in which the relation between the excursion set 



model and the conditional mass function are accounted for 
more carefully, are almost always even more accurate. 

2.4 Dependence on local density 

Following Mo & White (1996), knowledge of the conditional 
mass function allows one to estimate how the distribution of 
dark matter haloes today depends on the average density in 
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Figure 4. Conditional mass functions in the SCDM simulations (symbols with error bars). The two sets of curves in each panel show 
parent haloes with masses in the range 1 < M/M, < 2 (upper curves) and 16 < M/M» < 32 (lower curves); the upper curves have 
been shifted upwards by a factor of ten for clarity. Dotted curves show the spherical collapse prediction (Lacey & Cole 1993), dashed 
curves show the distribution one gets by rescaling the unconditional mass function (equation ^), and solid curves show our analytic 
approximation to the random walk with moving-barrier simulations (equation M). 



which the haloes are. In essence, they argued that a dense 
region should be thought of as an object which will collapse 
and form a virialized halo in the near future. This means 
that the haloes in it today can be thought of as 'progenitor 
subhaloes' viewed at 'low redshift'. In contrast, it will be 
a much longer time before an underdense region collapses 
(if it collapses at all) , so the haloes within it today are like 
the progenitor haloes viewed at high redshift. In hierarchical 
models, massive haloes form later, and less massive haloes 



form earlier. The discussion above means that one expects 
haloes in dense regions to be more massive, on average, than 
in underdense regions. The precise dependence of halo mass 
on local density depends on the precise relation between the 
local density today and the effective 'redshift'. Mo & White 
used the spherical collapse model to provide this relation. 
They provided a simple fitting formula to this relation in 
an Einstein-deSitter universe. Lemson & Kauffmann (1999) 
and Sheth & Lemson (1999a) showed that this provided a 
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Figure 5. Mass functions as a function of local density in the 
SCDM simulations (symbols with error bars), plotted as a func- 
tion of relative mass, so as to resemble the conditional mass 
functions presented earlier. Dotted curves show the spherical col- 
lapse prediction, dashed curves show the distribution one gets by 
rescaling the unconditional mass function (equation ^ , and solid 
curves show our analytic approximation to the random walk with 
moving-barrier simulations (equation The curves have been 
offset upwards by a factor of ten and a hundred, in the case of the 
middle and topmost curves, respectively. The upper most curves 
show the densest cells. 



reasonably good description of how, in their simulations, the 
density of haloes depended on local density. We have checked 
that the following simple modification to their formula is 
reasonably accurate for all cosmologies of interest: 



5o{S, 20 ) 



1.68647 



1.68647 



1.35 



(1 + 5)2/3 



1.12431 



+ 



0.78785 



(1+ ,5)1/2 (1 + 5)0.58661 



, (8) 



where M/pV = (1 + 5) is the nonlinear density of a re- 
gion containing mass M within the volume V at zq, and 
5sc{zo) denotes the critical density for spherical collapse 
at zo- The number density of m-haloes in regions of non- 
linear density S is got from the conditional mass function 
jV[m,SBc{zi)\M,5o{S, zo)] of equation 

The previous section showed that, in fact, the condi- 
tional mass functions are better fit by ellipsoidal collapse 
based curves. Therefore, one might reasonably expect that 
the same will be true for n{m\S). To emphasize how similar 
n{m\5) is to the conditional mass function, we have chosen 
to do the following. We have divided the simulation volume 
up into into cubes, each 10Mpc//i on a side. We then divided 
the cubes into three classes: the densest, and least dense ten 
percent of the cells, and the ten percent around the median 
density. Figs. |^ and ^ show n{m\5) for the haloes in the 
cells, plotted as a fraction of the mass of a cell. Symbols 



Figure 6. Same as the previous figure, but for ACDM. 

with error bars show the measurements in the cosmological 
simulations, histograms show the associated random walk 
distributions, solid curves show our analytic approximation 
to the random walks, dashed curves show the result one 
gets by rescaling the unconditional mass function, and dot- 
ted curves are for the spherical collapse, constant barrier 
model. The curves for the three types of cells have been off- 
set from each other for clarity; the lowest density cells are 
the lowest curves, and they have not been shifted, average 
density cells have been shifted upwards by a factor of ten, 
and the densest cells have been shifted upwards by a factor 
of a hundred. In all cases, the random walk model provides 
the best fit to the simulation data. 

Notice how similar these curves appear to the condi- 
tional mass functions presented earlier. Because most of 
haloes in the densest cells are a significant fraction of the 
total mass in the cell, the mass function in dense cells looks 
very like the low redshift conditional mass functions. In con- 
trast, the mass function in underdense cells looks much more 
like the high redshift conditional mass functions. This is pre- 
cisely what the model predicts. 

Figs. 1^ and ^ show what this trend means for the ac- 
tual number density of haloes in dense and less dense re- 
gions. The various symbols and curves are the same as in 
the previous two figures; the only difference is that now 
the x-axes have been multiplied by the total mass in the 
cell to show physical, rather than relative masses. Clearly, 
less dense cells have essentially no massive haloes. In addi- 
tion, the ratio of massive to less massive haloes is higher 
in denser cells. In particular, note that the density of less 
massive haloes in dense regions is actually smaller than the 
density of less massive haloes in underdense regions. One 
might have thought that dense regions simply have more 
haloes on average. For example, one might have thought 
that n{m\5) = (1 + 5)n{m). Our figures show that this is 
wrong, but that a good estimate of n{m\5) can, nevertheless, 
be computed analytically. We conclude this section with the 
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Figure 7. Same as Fig. but now presented as a function of 
physical, ratlier tiian relative, mass. Massive haloes occupy the 
densest cells (upper most curves). 



observation that our moving barrier based formulae provide 
a more accurate fit to the simulations than does the spheri- 
cal coUapse based constant barrier model. 



2.5 Rescaling the conditional mass function 

In the excursion set model with a constant barrier height, the 
unconditional mass function, when expressed as a function 
oi V = Sc/a'^i'm), is expected to be a universal function 
which is independent of redshift, cosmology or initial power 
spectrum. In addition, if the conditional mass function is 
expressed in units of [Sci — dco)'^ / {s — S), then it is expected 
to have the same shape as the unconditional mass function. 

The previous section we noted that, although the un- 
conditional mass function is a universal function of i/, this 
function is not the one predicted by the constant barrier 
model. We argued that if we interpret the unconditional 
mass function as coming from a moving barrier, then we no 
longer expect the conditional mass function to be a univer- 
sal function of u. Figs. ^ and ^ show this explicitly: they 
show the result of applying this rescaling to the conditional 
mass functions in the SCDM and the ACDM simulations we 
presented earlier. 

The four panels in each figure show the conditional mass 
function at each of the four redshifts we have been studying 
so far. The symbols show the result of rescaling the con- 
ditional mass functions in the simulations for parent haloes 
with mass in the range 1-2 (filled circles) and 8-32 M« (open 
triangles) at z = 0. Notice that the symbols in each panel do 
not overlap exactly — at fixed z, the conditional mass func- 
tions for different parent haloes do not rescale exactly. In ad- 
dition, the band traced out by the symbols at z = 4 is quite 
different from the band traced out at z = 0.5; the mass func- 
tions at different output times do not rescale either. Both 
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Figure 8. Same as the previous figure, but for ACDM. 

these findings illustrate our main point: the conditional mass 
function is not a universal function of the scaling variable u. 

The dotted curves which are the same in all the panels 
show the predictions of the constant barrier model; they do 
not provide a good fit at any time for any mass range. The 
dashed curves which are also the same in all the panels show 
the result of assuming that, upon rescaling, the conditional 
mass function will have the same shape as the unconditional 
mass function; although they provide a good fit at large z, 
they are increasingly in error at small z. This is true both 
for the SCDM and the ACDM models. 

The solid lines in the various panels show the predic- 
tions of our moving barrier model. In this case, the predic- 
tions depend both on the parent mass range, and on the red- 
shift at which the progenitors are identified: we have chosen 
to show the predictions for the 1-2M* haloes only. Whereas 
the model is in reasonable agreement with the simulations 
at large z, it has the wrong shape an small z. This is not 
terribly surprising, because our formula for the first crossing 
distribution of the moving barrier model was only supposed 
to work at large z, but it is unfortunate that the disagree- 
ment at low 2 is so bad! One might have thought that the 
actual first crossing distribution may be in good agreement 
with the GIF simulations, and that it is only the analytic ap- 
proximation which is in error at small z. Unfortunately, this 
is not so. The histograms in Fig. ^ show the result of sim- 
ulating an ensemble of random walks to construct the con- 
ditional mass functions; although they are in slightly better 
agreement with the simulations, they are still quite different. 

The discrepancy between the simulations and our mov- 
ing barrier model predictions are most pronounced when 
the subclump mass is m/M > 1/2. This suggests that our 
model is unable to describe the histories of clumps at small 
lookback times. At small lookback times, one might worry 
that the spherical overdensity groupfinder we use to iden- 
tify the subclumps in the simulations might find a different 
set of objects than a friends-of-friends algorithm. Plots of 
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Figure 9. Rescaled conditional mass functions in the SCDM model. Panels show the different redshift bins we studied earlier. Symbols 
in each panel show the different mass ranges we considered. Dotted line shows the constant barrier prediction (in these variables it is 
the same as the unconditional mass function), dashed line shows the result of rescaling the actual unconditional mass function, and solid 
curves show the result of rescaling the moving barrier predictions. 



the rescaled conditional mass function constructed using a 
friends-of-friends algorithm look very similar to the spheri- 
cal overdensity results presented above — the discrepancy be- 
tween model and simulations is independent of the choice of 
groupfinder. 

In addition, because the theory assumes that mass is 
conserved — all the mass of a subclump becomes part of the 
final halo — whereas this is not true in the simulations: some 
of the particles which make up the final object may have 
come from a subclump which merged along with most of 



its particles into a different object. This means that there 
is some freedom associated with how we decide if a clump 
at an early time should be counted as a subclump of a halo 
at the final time. We have tried two schemes for idenfifying 
subclumps: a progenitor is a clump which donates at least 
half its mass to the final object, or which donates at least one 
particle to the final object. Once we have made this decision, 
we must also decide what we wish to count as the mass of the 
parent object: two natural choices are the mass at the final 
time (which, by definition is fixed for all earlier redshifts) , or 
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Figure 10. Same as the previous figure, but for ACDM. 



the mass which is got by summing up the masses of all the 
progenitor subclumps (which may depend on redshift). The 
figures above are for the case in which a subclump is any 
clump which donates at least one particle to the final ob- 
ject, and the parent mass is defined as the mass at the final 
time (so it is independent of redshift). Although the actual 
conditional mass functions depend slightly on which com- 
bination of the above choices we make, the generic results 
shown above are independent of this choice. 

Before moving on, we think it worth noting that the dis- 
crepancies between the SCDM simulations and the dotted or 
dashed curves are qualitatively similar to the discrepancies 



in the ACDM case. This suggests that one should be able 
to find a model which can account for these discrepancies in 
a way which is independent of power-spectrum, redshift or 
cosmology. Our moving barrier model is just not up to the 
task. The next section studies why. 



3 AN EXTENSION 

In the ellipsoidal collapse model envisaged by Bond & Myers 
(1996) and implemented by Sheth, Mo & Tormen (2001), 
the collapse of a patch is determined by the surrounding 
shear field. In a Gaussian random field, the field around a 
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Figure 11. Distribution of overdensities and scales at which the 
six-dimensional random walks crossed the ellipsoidal collapse bar- 
rier. Solid curve shows the approximation used by Sheth, Mo 
& Tormen (2001) in their, considerably simpler, one-dimensional 
random walks. 



patch may differ from patch to patch. Appendix |B| provides 
a simple prescription for choosing a set of patches which 
have the correct ensemble averaged properties — in essence, 
this requires studying the first crossing distribution of six- 
dimensional random walks. 

Because a six-dimensional walk is computationally ex- 
pensive, rather than choosing the distribution of initial 
patches from this exact statistical distribution, Sheth, Mo & 
Tormen suggested it should be a good approximation to use 
an appropriately chosen mean value, and neglect the scat- 
ter around this value. This allowed them to reduce what is 
a six-dimensional random walk to a one-dimensional walk. 
It is this one-dimensional walk which we have considered 
so far. One might worry, however, that the discrepancy be- 
tween model predictions and simulation results we found in 
the previous section may actually be due to our neglect of 
the scatter around the average value. The main purpose of 
this section is to study this possibility. 

Fig. ^ shows the result of generating an ensemble of 
four thousand six-dimensional random walks associated with 
Gaussian random fields as described in the Appendix. The 
walks are stopped when they cross the barrier associated 
with the ellipsoidal collapse model of Bond & Myers (1996). 
The crosses show the values of 5 and a at which the six- 
dimensional walks crossed the ellipsoidal collapse barrier 
5crit(e,p): we actually used the simple fit, equation 3 in 
Sheth, Mo & Tormen (2001), to the critical density required 
for collapse 5crit(e,p). The solid curve shows the approxima- 
tion used by Sheth, Mo & Tormen; it provides a rather good 
description of how 5crit increases with increasing a. The Ap- 
pendix describes the reason for this in more detail. For now, 
we simply note that because the solid curve provides a rea- 
sonably good description of the crosses, the first crossing 
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Figure 12. Scaled unconditional (jagged curves) and conditional 
(filled circles) mass functions associated with six-dimensional ran- 
dom walks which cross the ellipsoidal collapse barrier. Dashed 
line shows the unconditional mass function which fits the mass 
function of bound objects in simulations of clustering well. As dis- 
cussed in the previous section, the dashed curve is well described 
by the distribution of first crossings by one-dimensional random 
walks, of a barrier which is associated with ellipsoidal collapse. 
Dotted line shows the corresponding spherical collapse prediction. 



distributions of the six-dimensional walks considered here 
are unlikely to be very different from the first crossing dis- 
tributions associated with the (considerably simpler) one- 
dimensional walks studied in the previous sections of this 
paper. 

The two jagged curves in Fig. ^ show this explicitly. 
They show the first crossing distributions associated with 
the six-dimensional walks which cross the z = and z — 0.5 
six-dimensional ellipsoidal collapse barriers. They have been 
rescaled similarly to how the unconditional mass functions 
in simulations rescale: v = {Sbc/o')^ . After rescaling, the two 
curves appear similar, as they should; note that they are 
reasonably like the dashed curve, which shows equation (^, 
and they are rather different from the dotted curve which 
shows the spherical collapse prediction. 

Recall that the dashed curve is very similar to the mass 
function one gets by simulating one-dimensional random 
walks. Because the jagged curves are in quite good agree- 
ment with the dashed curve, the first crossing distributions 
associated with the six-dimensional walks are actually rather 
similar to those associated with the (considerably simpler) 
one-dimensional walks studied in the previous sections. This 
agrees with what our conclusions from Fig. 0. The Ap- 
pendix describes the reason for this in more detail. 

Ifaving shown that the unconditional mass functions as- 
sociated with the six-dimensional random walks are in rea- 
sonable with numerical simulations, we now turn to the con- 
ditional mass functions — the test which our one-dimensional 
random walks failed. The solid symbols with error bars 
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show the conditional mass functions associated with the six- 
dimensional walks, expressed in the scaled units of the pre- 
vious section: u = 0.707 S^^ 0.5^ / (s — S), where s denotes the 
value of at which the z = 0.5 six-dimensional boundary 
was crossed, and S is the scale on which the lower z = 
six-dimensional boundary was crossed. Note the excess of 
points at large values of v; these conditional mass functions 
are quite different from the conditional mass functions in 
simulations. Indeed, the discrepancy between the excursion 
set predictions and the simulations of hierarchical clustering 
has got worse! 

There are at least two reasons why the excursion set 
approach with one- dimensional random walks may fail at 
small lookback times. The first is the excursion set neglect 
of correlations between scales (Peacock & Heavens 1990; 
Bond et al. 1991). At large lookback times most subclumps 
are a small fraction of the mass of the parent halo, so the 
smoothing scale associated with the subclumps is sufficiently 
different from that of the parent that the neglect of correla- 
tions between the two scales is probably justified. At smaller 
lookback times the parent and subclump scales are not so 
well separated, so the neglect of correlations is a more likely 
to be a bad approximation. This is one possible reason for 
the agreement at large lookback times despite the discrep- 
ancy at low redshift. The second possibility is that the one- 
dimensional parametrization of ellipsoidal collapse outlined 
by Sheth, Mo & Tormen (2001) is too simple. The resuhs 
of this section suggest that, in fact, it is the first possibility 
which is the cause of the discrepancy. 



4 DISCUSSION 

Sheth, Mo & Tormen (2001) argued that a simple modifi- 
cation to the original excursion set approach was enough to 
improve agreement between the predictions of the approach 
and numerical simulations. The modification they suggested 
was to the value of the linearly extrapolated critical over- 
density 5c associated with the collapse of an object. The 
spherical collapse model assumes that this value is indepen- 
dent of the mass of the collapsed object, whereas ellipsoidal 
collapse makes 5c depend on m. In the context of the ex- 
cursion set approach, this corresponds to studying the first 
crossing statistics of a set of moving, rather than constant 
barriers. We also argued that a moving barrier also provides 
a simple way in which the excursion set appraoch can be ex- 
tended to apply to models in which the initial dark matter 
distribution is not completely cold. 

Because moving barrier models are so useful, we pro- 
vided analytic approximations for the required first crossing 
distributions (equations ^ and |^ , and showed that they were 
reasonably accurate. Although our formulae for the condi- 
tional mass functions (solid curves in Figs. |^-|^) are slightly 
different from, and usually more accurate than, those one 
obtains by a simple rescaling of the unconditional mass func- 
tion (dashed curves in the same figures), this simple rescal- 
ing of the unconditional mass function is still more accurate 
than what one gets if one rescales the constant barrier for- 
mulae (dotted curves) . We showed that the predicted uncon- 
ditional and conditional mass functions were in reasonably 
good agreement with results from numerical cosmological 
simulations (Figs. pHsl), and we also provided a simple effi- 



cient algorithm which allows one to generating masses which 
have the correct universal mass function (Section 2.z). 



However, we showed that neither the constant nor the 
moving barrier models were able to describe the simulation 
results at small lookback times (Figs, ^and ^o|). This means 
that our results for the conditional mass functions cannot 
be used to generate realizations of the forest of merger his- 
tory trees. We argued that this discrepancy was most likely 
due to the excursion set approach's neglect of correlations 
between scales (e.g. Peacock & Heavens 1990; Bond et al. 
1991; Monaco 1997b). While this neglect is a bad approx- 
imation at small lookback times, it is reasonably accurate 
at large lookback times. This is why the excursion set ap- 
proach is able to provide a reasonably good description of 
clustering at high redshift, even though it is inaccurate at 
small redshifts. 

Before concluding, we will consider how some of our re- 
sults are related to other work in the literature. Recently, 
Jenkins et al. (2001) showed that, although the mass func- 
tions in their simulations scaled in accordance with the ex- 
cursion set prediction, our equation (^) slightly overesti- 
mated the unconditional mass functions in their simulations. 
We thought it would be interesting to show the various ap- 
proximations to the mass function which we presented in this 
paper on one plot. The dot-dashed curve in Fig. ^ with 
cutoffs at low and high masses shows the fitting function 
Jenkins et al. proposed, which fits their simulations well, the 
dashed curve shows the fitting formula of equation (^ with 
a = 0.707 (following Sheth & Tormen 1999), the histogram 
shows the distribution one gets by simulating random walks 
with the ellipsoidal collapse moving barrier (equation |l], fol- 
lowing Sheth, Mo & Tormen 2001), the solid curve shows 
our approximate formula for this first crossing distribution 
(equation ^ , and the dotted curve shows the spherical col- 
lapse, constant barrier prediction (Press & Schechter 1974; 
Bond et al. 1991). Simulations currently available do not 
probe the regime where u < 0.3 or so (the Jeans mass is at 
about 1/ ~ 0.03). 

The upper set of curves show the residuals between 
our formulae and the one provided by Jenkins et al., in the 
regime to the right of the low-mass cut-off (marked by an ar- 
rowhead). In addition to the previously mentioned formulae, 
we have included a new dashed curve which shows the result 
of changing a in equation (^ from 0.707 to 0.75. This simple 
change appears to be all that is necessary to reduce the dis- 
crepancy between it and the simulations substantially. Our 
formula differs dramatically from the one Jenkins et al. pro- 
pose at small masses. We hope that simulations in the near 
future will be able to address which low mass behaviour is 
correct. 

To illustrate that our formulae really do work for a large 
class of barrier shapes. Fig. |lj shows the first crossing dis- 
tribution associated with the barrier discussed by Monaco 
(1997a,b): 



B{S) 



1.82/(5sc - 0.69 y/s/sl 



where, as throughout this paper, S* = 5sc- The height of 
this barrier decreases with S, so it really is quite different 
from ours (Sheth, Mo & Tormen 2001 discuss the physical 
reason why). (We chose not to present results for the bar- 
rier shape studied by Del Popolo & Gambera (1998) because 
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Figure 13. Comparison of the various mass functions described 
in this paper, with the one which fits the Hubble volume simu- 
lations presented in Jenkins et al (2001). Upper curves show the 
residuals between our analytic formulae and the Jenkins et al. 
fitting formula. 
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their shape is not so different from ours, whereas Monaco's 
really is quite different. The fact that this barrier decreases 
with S means that all walks are guaranteed to cross it, and 
that there is no fragmentation associated with this barrier 
shape. (This is in contrast to barriers whose height increases 
sufficiently strongly with 5"; for the Sheth, Mo & Tormen 
2001 barrier shape studied in the main text, unbound mass 
and fragmentation are features which are formally possible 
but rare in practice.) The histogram shows the numerical 
Monte-Carlo first crossing distribution, and the two solid 
curves shows our analytic approximation, computed by in- 
serting this barrier shape into our equation (U|). The curve 
which provides a slightly worse fit to the histogram shows 
the result of using the first five terms in the series (as we did 
for the other figures in this paper); the other curve shows 
what happens if we use the first ten terms instead. Just for 
comparison, the dotted and dashed curves show the spher- 
ical collapse prediction, and the one which actually fits the 
cosmological simulations (equation ^ . The figure shows that 
we are able to describe the first crossing distribution of this 
barrier shape well. This means that one could, in principle, 
use our formula, with Monaco's barrier, to study the con- 
ditional mass functions associated with his parametrization 
of nonlinear collapse — we have not pursued this further, al- 
though comparison of this first crossing distribution with 
the z — 0.5 panels in Figs. ^ and ^ suggest that this might 
be a useful exercise. 

In summary, we have provided a formula which de- 
scribes the first crossing distribution of independent uncor- 
related Brownian motion random walks, for a wide class of 
moving barriers. This formula can be used to provide simple 
accurate predictions for a number of statistical quantities as- 
sociated with the formation and clustering of dark matter 
haloes, all within the same formalism. 
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Figure 14. The unconditional mass function associated with 
a model of collapse which was presented in Monaco (1997a, b). 
Whereas the height of the barrier studied in our paper increases 
with distance, the one proposed by Monaco decreases with dis- 
tance. Despite the dramatic differences between the two barriers, 
inserting this shape into our equation (^ provides a good fit (solid 
curve) to the exact result (histogram). 
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APPENDIX A: THE LINEAR BARRIER 

Suppose that the barrier shape increases linearly with in- 
creasing variance S = a^: 

B{S,z)=5,{z)(^l + S/S4z)y (Al) 

where Sc{z) = (5co(l + z) and S*(z) = Sto{'^ + z)^ if SI = 
1. These scalings with z are what is required by the self- 
similarity of Brownian motion. In the Bond et al. (1991) 
formulation of the constant barrier problem, it is customary 
to express the S axis in the units it had at some fiducial 
time, say z = 0, and to study the successive crossings of 
barriers having different values of z. In effect, the constant 
height Press-Schechter barrier has S*o = cxd, so in that case 
only the scaling of the i/-axis was apparent. 

This linear barrier shape is motivated by the observa- 
tion that the GIF simulations have fewer low mass haloes 
relative to high mass ones, as compared to what is predicted 
by the constant barrier model. This means that, at least for 
some range of S, the moving barrier must have a positive 
slope, since this would make it relatively easier to cross at 
small S (large mass) than at large S (small mass), as com- 
pared to a barrier of constant height. 

An additional reason for considering this linear barrier 
is the following. Recent work (Bode, Ostriker & Turok 2001) 
suggests that the halo mass function in which the dark mat- 
ter is warm initially has even fewer low mass halos than cold 
dark matter based ellipsoidal collapse models predict. In the 
context of the approach outlined by Sheth, Mo & Tormen 
(2001), the physical reason for this is relatively simple: low 
mass haloes do not form because they are hotter initially, so 
a larger overdensity is required to hold them together against 
the thermal pressure which prevents collapse, or against the 
stronger shearing from the velocity field. This suggests that 
the critical density required for collapse by the present time, 
Sccim), should increase even more strongly with decreasing 
m than it does when the dark matter is cold. The warm 
dark matter model is not particularly well motivated, and 
its free parameters have not yet been fixed, it seems prema- 
ture to provide a detailed S^°^{m) relation at the present 
time. For this reason, the linear barrier considered in this ap- 
pendix should be thought of as an example of what happes 



when the barrier height increases even more steeply with 
decreasing mass than it does in the Sheth, Mo & Tormen 
cold dark matter models. 

In the constant barrier model, all random walks were 
guaranteed to cross the barrier. This is because the rms 
height of random walks at S is proportional to y/S, so at 
sufficiently large 5", all walks will have crossed the constant 
barrier. As a result, the associated first crossing distribution 
is normalized to unity: since each random walk is associated 
with a volume element in the initial Lagrangian space, this is 
usually interpretted as meaning that, in the constant barrier 
model, all the mass in the universe is bound up in collapsed 
objects of some mass, however small. In contrast, the lin- 
ear barrier ( |A1| ) increases to arbitrarily high values at high 
5*. Because the rms height of the random walk grows more 
slowly than the rate at which the barrier height increases, 
there is no guarantee that all random walk trajectories will 
intercept this barrier. Indeed, only a fraction e~'^''<: '^''''^*'^'' 
of them will. So, in the linear barrier model, not all initial 
volume elements are associated with bound haloes. Since not 
all particles in numerical cosmological simulations are asso- 
ciated with bound haloes anyway (the fraction of unbound 
mass is typically on the order of ~10%, though how much of 
this is a consequence of limited resolution in the simulations 
is uncertain), this feature of the moving barrier model may 
or may not be a good thing. In any case, this is one qual- 
itative difference between the moving boundary model and 
a model with constant barrier height. (Readers who dislike 
this feature of the linear model are invited to patch a con- 
stant barrier of height 5co at small 5 to a linearly increasing 
one at intermediate S to another constant one (but now at 
at a greater height, of course) at a large S of their choice, 
and to compute the associated mass function! This can be 
done relatively easily using the results presented below.) 

In addition to the question of normalization, barriers 
whose heights increase with 5* more rapidly than 5°'^ will 
have mass functions in which the low mass end is depleted 
relative to the constant barrier case. In the linear barrier 
model studied below, the mass function has an exponen- 
tial cutoff at both low and high masses. Because the barrier 
shape associated with ellipsoidal collapse grows only as S"'® 
both these effects are extremely weak. Nevertheless, it is 
worth bearing these features in mind as bigger and better 
simulations become available. 

The second important qualitative difference between 
a moving barrier model and the original constant barrier 
model is the following: whereas the y-intercept B{0,z) in- 
creases as z increases, the slope decreases as {l + z)~^. This 
means that it is possible for barriers to intersect at finite 
values of S. For example, two linear barriers B{S, zq) and 
B{S,zi), with zi > zq, will intersect at that critical value of 
S at which B{S,zo) =B{S,zi): Soi/S,o = (l + zo){l + z^). 
This means that all trajectories which first cross -6(20) at 
5* > Soi must have crossed -8(^1) at a smaller value of 5*. The 
logic of Lacey & Cole (1993) then says that all haloes at zo 
that are less massive than the associated critical mass Moi 
were formed by fragmentation of a halo which, at zi > zq, 
was more massive. 

Again, this may or may not be a good thing. Simula- 
tions show that ~20% of the total mass ever associated with 
progenitors of a halo does not find its way to the parent halo 
(Tormen 1998). Presumably this reffects the fact that while 
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Figure Al. The conditional mass functions associated with the linear barrier. Smooth solid curves show the analytic linear barrier 
prediction, dashed curves show the conditional constant barrier distribution for comparison, and histograms show the result of our 
numerical Monte— Carlo calculation. 



small haloes may fragment, the fragmented mass is not a 
very large fraction of the total mass. Simulations also show 
that some subclumps pass through the virial radius of a 
given parent halo a number of times, each time depositing 
some fraction of their mass, before they finally become part 
of the parent halo. It may be that a moving barrier model is 
able to incorporate and perhaps even quantify these effects. 

The first crossing distribution of trajectories that do 
cross the linear barrier is Inverse Gaussian: 



f{S,z) dS = 



B{Q,z) 



exp 



2S 



AS_ 

S ' 



(A2) 



(Whitmore 1978; also see Sheth 1998). The associated un- 
conditional mass function is got by inserting this in equa- 
tion (^. Whereas the mass function in the constant barrier 
case has an exponential cutoff at the high mass end, the 
mass function in this linear barrier model has exponential 
cutoffs at both low and high masses. Recall this arises from 
the fact that barrier height increases so rapidly at large 5*. 
Therefore, the mass function associated with this linear bar- 
rier has fewer small mass objects than the constant barrier 
predicts, in qualitative agreement with the GIF simulations. 
(The agreement is certainly not quantitative, but the linear 
barrier is only used here for illustrative purposes.) 

Since a linear barrier is linear whatever the origin of the 
coordinate system, the solution to the two barrier problem 
is also Inverse Gaussian. That is, given zi > zo and given 
that the first crossing of B{zo) occured at So < 5*01, the 
probability that the first crossing of B{zi) occurs in the 



range dSi about 5*1 is given by equation (AS), with the 
substitutions S {Si — So), and B Bio, where 

B-w{Si - So) = B{Si,zj) - B{So,zo). 



As was the case for the unconditional mass function, the 
height of this barrier diverges as 5i — 5*0 — > oo, so not all 
trajectories intersect it. Again, it seems reasonable to as- 
sociate the fraction that do not with the fraction of the 
parent halo mass that is not associated with bound sub- 
clumps. (Those readers who computed the mass functions 
associated with the patched constant-linear-constant barri- 
ers may disregard the previous two sentences, provided they 
first compute the associated two-patched-barrier problem!) 

show this condi- 



The smooth solid curves in Fig. Al 



tional distribution for a few representative choices of the 
parent halo mass M. The parents were assumed to have 
formed ai zq = 0, and the progenitor distributions are shown 
at an earlier redshift zi . The underlying power spectrum was 
chosen to be the same as the GIF SCDM power spectrum. 
It is also straightforward to solve this problem numerically: 
for the values of M shown in the figures, the histograms 
show the distribution generated by simulating the crossings 
of the higher linear barrier by 10'' random Brownian motion 
trajectories that started at the initial positions B{S, z = 0), 
where S{Ad) is given by the GIF power spectrum. This fig- 
ure has been included mainly to show that our Monte-Carlo 
code works, because we will use the code in the next section. 
We have also verified that the numerical code gives the cor- 
rect conditional and unconditional mass functions when the 
barrier heights are constant. 

For completeness, we also show the bias relations asso- 
ciated with this model. Following Mo & White (1996) (also 
see Sheth & Tormen 1999), the mean Lagrangian bias be- 
tween haloes and mass is 



f{Sl,Zl\So, Zq) 

/(5i,zi) 



(A3) 
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The limit Mq ^ Mi and \So\ <C 5i is sometimes called the 
peak-background split. In this limit 

(5h(l|0) — > (5o, where fi = dc{zi)/Si 

for the linear barrier model. Massive haloes have small values 
of S, so they have large values of v. Small haloes have v ^ 0. 
Thus, in this limit, less massive haloes are unbiased relative 

to the mass, whereas massive haloes arc positively biased. 
For comparison, the corresponding limit for the constant 
barrier is 

Thus, the predictions of the linear and constant barriers are 
similar for massive haloes, but they differ for less massive 
ones. In particular, less massive haloes in Lagrangian space 
arc anti-biased in the constant barrier model, whereas they 
are unbiased in the linear barrier model. 

The halo-to-mass bias in the evolved Eulerian space for 
the constant barrier can be computed by expanding 

5^{l\0) = {l + 5) [l+Stim]-! (A4) 

to lowest order in 5 (Mo & White 1996). In this limit, 5 m So, 
so (5j^(l|0) = (1 -I- [vi — l]/5i)S for the constant barrier. 
The same logic gives (5j^(l|0) = (1 + ui/5i)5 for the linear 
barrier. Again, the predictions of the constant and linear 
barriers agree for massive haloes, but are different for less 
massive ones; the haloes in the linear barrier model are more 
positively biased. This is encouraging, because, as mentioned 
in the introduction, this is in qualitative agreement with the 
trend seen in numerical simulations. 

The rate of increase of the ellipsoidal collapse barrier 
is shallower than for the linear barrier discussed here. Since 
it is intermediate between the linear barrier and the con- 
stant spherical collapse barrier, we might reasonably expect 
the laxge scale bias factor of less massive haloes in the el- 
lipsoidal collapse model to be greater than that associated 
with the spherical collapse model. Since the ellipsoidal col- 
lapse barrier rises less steeply than linear, the modified large 
scale bias will be somewhat less than the linear barrier pre- 
diction. Fig. 4 of Sheth, Mo & Tormen (2001) shows that 
this is, indeed, the case. 

We argued at the start of this section that warm dark 
matter models can be parametrized by making the criti- 
cal density for collapse depend more strongly on mass than 
when the dark matter is cold. In this respect, a comparison of 
the linear barrier formulae given here with the results for the 
barrier written down by Sheth, Mo & Tormen (2001) shows 
why, generically, warm dark matter models are expected to 
have fewer low mass haloes and, consequently, different bias 
relations, paxticularly at the low mass end, than cold daxk 
matter models. This is in qualitative agreement with the 
numerical simulations of Bode, Ostriker & Turok (2000). 



APPENDIX B: DISTRIBUTION OF DENSITY 

AND ANGULAR MOMENTUM OF A PATCH 
IN A GAUSSIAN RANDOM FIELD 

Let dij = Vijtl), where 4> is the initial potential, denote the 
various components of the deformation tensor D (here 1 < 
i < 3 and similarly for j). Following, e.g. Bardeen et al. 
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(1986), these components are 

dii = (-2/1 -3t/2/Vl5-j/3/V5)/3 

d22 = {-yi + 2j/3/\/5)/3 

d33 = (-t/i+3j/2/\/l5-t/3/\/5)/3 

di2 = yt/VTs 

d23 = y5/VT5 

di3 = ye/VTE (Bl) 

where the yiS are independent Gaussian variates with zero 

mean and variance . Poisson's equation says that the trace 
of this matrix, Tr(I?) equals the overdensity 5. Thus, 

5 = dii +d22 + d33 = -yi; (B2) 

the final expression shows explicitly that 5 is a Gaussian 
random variate. 

The eigenvalues of this matrix are the roots of the char- 
acteristic equation 

P(A) = I)et[V - AX] = -oo - oiA - 02A^ - A^ (B3) 

where T is the identity matrix, and the o^s are various com- 
binations of the dijS got by expanding the expression above 
and ordering by powers of A. Since I? is a 3 x 3 real sym- 
metric matrix, P(A) is a cubic with three real roots which 
satisfy 

Al -I- A2 -I- A3 = —a2 
A1A2 + A2A3 -|- A1A3 = ai 

A1A2A3 = -00. (B4) 

Because rotations leave the trace unchanged, the overden- 
sity 5 is the sum of the three eigenvalues, so — 02 = S = —yi- 
In addition, the square of the angular momentum is propor- 
tional to 

2 ^ (Ai-A2)^ (A2-A3)' (Ai-A3)^ 
2 2 2 

= 0,2 — 3ai (B5) 

(e.g. Heavens & Peacock 1988; Catelan & Theuns 1998). 
This shows that to get we don't need to solve the cubic — 
we just need to read off the appropriate coefficients of the 
characteristic equation. Thus, we find that 

= + 3(d?2 + + 1^23 — ^11^22 — diidss — ^22^33) 
= (j/2 +2/3 +1/4 +2/5 +J/6)/5 (B6) 

Although the first line suggests that is coupled to S, the 
final expression shows that it is not. In particular, the ex- 
pressions above show that S is distributed as a Gaussian, 
and is an independent variate drawn from a Chi-square 
distribution with five degrees of freedom, xii'^)- The fact 
that the overdensity and the square of the angular momen- 
tum are independent does not seem to have been noticed 
before. A XsC"") distribution is rather similar in shape to a 
Lognormal, so this provides a simple way to see why spins 
of peaks in Gaussian random fields are also approximately 
lognormal (Heavens & Peacock 1988). 

The results above can be used to generalize the excur- 
sion set algorithm studied in the main text. Set n = and 
yi{n) = for 1 < i < 6. Thereafter, at each step labeled by 
n, choose six, rather than one, independent Gaussian ran- 
dom variates gi, each with variance s. For each i = 1, 6 set 
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Figure Bl. Comparison of two prescriptions for collapse, both of 
wliicli give good agreement with the unconditional mass function. 
Crosses show the actual values of collapse densities and scales, 
and solid curves show the approximation assumed by Sheth, Mo 
& Tormen (2001). 



yi{n) — yi{n — I) + Qi, and use these to compute 5 and r^. 
These can be used to give the values of the overdensity and 
the angular momentum for the scale on which the variance 
is oc ns (recall that n is the number of steps taken). 
Now check to see if 5 = —yi{n) exceeds a critical value, 
say ScTit{o',r^)- If it does, the six-dimensional walk stops at 
this scale. If not, the walk continues to smaller scales. Be- 
cause each of the giS is chosen independently of the values 
of the j/iS or of a, the walk takes independent steps in the 
six-dimensional space; it is in this sense that this algorithm 
generalizes the one-dimensional excursion set random walk 
studied in the main body of this paper. 

The analysis above shows clearly that the one- 
dimensional random walk approach of Sheth, Mo & Tormen 
(2001) corresponds to the following approximation. Replace 
the dependence of 5crit(o", r^) on the random variate by 
a dependence on its average value (r^) oc a^. This means 
that the critical density for collapse is a function of a alone, 
ScTit{cr). As a result, the random walk in six-dimensions can 
be reduced to a walk in one-dimension only, thereby greatly 
reducing the complexity of the problem. 

Recently just such a six-dimensional random walk al- 
gorithm has been used by Chiueh & Lee (2001), although 
they did not notice the considerable simplifications which 
follow from the algebra presented above. They simulated 
an ensemble of six-dimensional random walks, and set the 
parameters of the barrier to be crossed by requiring that 
the resulting first crossing distribution give the uncondi- 
tional mass function. In particular, they showed that S^it ~ 
1.5[1 + {2r^ /3f /0.15]°-'^^ provided a good fit to the required 
critical value of the overdensity. 

The algebra above allows one to see what such an ap- 
proach implies. To do this, suppose the barrier is 5crit = 



1.686(1 + r^). The dependence on means that if the par- 
ticle has walked to S = 1.686, it will still not have crossed, 
because is always certainly greater than zero. So, to 
cross, the particle needs to have some S > 1.686. How much 
greater? This depends on the typical value of r^. Because 
is drawn from a xi distribution, (r^) ~ a^. Now suppose 
that the xi distribution is very sharply peaked at its mean 
value (it quite well peaked, but taking the extreme case helps 
to see the argument). This means that the barrier shape is 
something like Sait = 1.686(1 + a^). The fact that a xi dis- 
tribution is not very sharply peaked at its most probable 
value simply means that sometimes when 5 = 1.686(1 +cr'^) 
the particle will still be less than (5crit, so the walk must go 
on. Of course, sometimes < cr^, and in this case the walk 
will stop even if 5 < 1.686(1 -I- a'^). So, this means that we 
can think of the dependence on as making the critical 
value of the boundary height, when expressed as a function 
of cr^ (the way Sheth, Mo & Tormen 2001 did) a little fuzzy. 
So, provided the xi distribution is not too broad, the con- 
siderably simpler one-dimensional random walk approach of 
Sheth, Mo & Tormen (2001) should be a reasonable approx- 
imation. 



Fig. Bl shows the result of doing this for two choices of 



the barrier shape, both of which produce first crossing dis- 
tributions which, when inserted into equation (^), give mass 
functions of bound objects which are similar to the one in 
simulations of hierarchical clustering. The upper panel was 
constructed using a barrier whose height increased linearly 
with r^, and the lower panel shows results for the barrier 
shape used by Chiueh & Lee (2001), which increases as 
r*. The crosses show the values of 5 and a at which each 
six-dimensional random walk crossed the barrier. The solid 
curve shows the approximation used by Sheth, Mo & Tormen 
(2001); it provides a reasonable description of the increase 
of 5crit with a. 

Determining the barrier shape by requiring agreement 
with the clustering simulations is unsatisfying, especially 
in view of the fact that the two different boundaries given 
above provide equally adequate approximations to the mass 
function. For this reason, the main text shows the result 
of combining the six-dimensional walk described above with 
the ellipsoidal collapse model of Bond & Myers (1996). This 
was relatively easy to do, because a simple fitting function 
for how the critical collapse boundary associated with this 
ellipsoidal collapse depends on the initial shear field has been 
given by Sheth, Mo & Tormen (their equation 3 and their 
Fig. 1). 
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